The PKSFC package: simulating simple models

Antoine Godin (Kingston University)

Agent-based and stock-flow consistent modelling: theory and applications - Paris - July 17

PKSFC package

Installing the dependent libraries and the package

You need to install all the required libraries This is for traditional libraries

install.packages("expm")
install.packages("igraph")

For non-conventional libraries, such as the one need to visualize Direct Acyclical Graphs (DAG), you need to do the following

source("http://bioconductor.org/biocLite.R")
biocLite("Rgraphviz")

Finally you can then download the PKSFC package from github and install it locally

install.packages("path/PKSFC_1.5.tar.gz", repos = NULL, 
                                 type="source")

Testing

Now we’re ready to load the package:

library(PKSFC)
## Loading required package: expm
## Loading required package: Matrix
## 
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
## 
##     expm
## Loading required package: igraph
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
  1. Load SIM (download SIM.txt from Github)
sim<-sfc.model("Models/SIM.txt",modelName="SIMplest model 
        from chapter 3 of Godley and Lavoie (2007)")
  1. Simulate the model
datasim<-simulate(sim)
  1. replicate figure 3.2 of page 73.
matplot(sim$time,datasim$baseline[,c("Yd","C")],type="l", 
    xlab="", ylab="",lwd=2,lty=c(2,3))
lines(sim$time,vector(length=length(sim$time))
    +datasim$baseline["2010","C"], lwd=2)
legend(x=1970,y=50,lwd=2, legend=c("Disposable Income",
    "Consumption","Steady State"),lty=c(2,3,1), bty="n")

How does it work?

Source code of SIM

#1. EQUATIONS
C = C
G = G
T = T
N = N
Yd = W*N - T
T = theta*W*N
C = alpha1*Yd + alpha2*H_h(-1)
H = H(-1) + G - T
H_h = H_h(-1) + Yd - C
Y = C + G
N = Y/W
#2. PARAMETERS
alpha1=0.6
alpha2=0.4
theta=0.2
#EXOGENOUS
G=20
W=1
#INITIAL VALUES
H=0
H_h=0
#3. Timeline
timeline 1945 2010

A few important points regarding the model source code:

Internal representation

print(sim)
## $name
## [1] "SIMplest model \n\t\tfrom chapter 3 of Godley and Lavoie (2007)"
## 
## $simulated
## [1] FALSE
## 
## $variables
##       name     initial value description    
##  [1,] "Yd"     NA            "1. EQUATIONS" 
##  [2,] "T"      NA            ""             
##  [3,] "C"      NA            ""             
##  [4,] "H"      "0"           ""             
##  [5,] "H_h"    "0"           ""             
##  [6,] "Y"      NA            ""             
##  [7,] "N"      NA            ""             
##  [8,] "alpha1" "0.6"         "2. PARAMETERS"
##  [9,] "alpha2" "0.4"         ""             
## [10,] "theta"  "0.2"         ""             
## [11,] "G"      "20"          "EXOGENOUS"    
## [12,] "W"      "1"           ""             
## 
## $endogenous
##      name  initial value lag description   
## [1,] "Yd"  NA            "0" "1. EQUATIONS"
## [2,] "T"   NA            "0" ""            
## [3,] "C"   NA            "0" ""            
## [4,] "H"   "0"           "1" ""            
## [5,] "H_h" "0"           "1" ""            
## [6,] "Y"   NA            "0" ""            
## [7,] "N"   NA            "0" ""            
## 
## $equations
##      endogenous value equation                 description   
## [1,] "Yd"             "W*N-T"                  "1. EQUATIONS"
## [2,] "T"              "theta*W*N"              ""            
## [3,] "C"              "alpha1*Yd+alpha2*H_h_1" ""            
## [4,] "H"              "H_1+G-T"                ""            
## [5,] "H_h"            "H_h_1+Yd-C"             ""            
## [6,] "Y"              "C+G"                    ""            
## [7,] "N"              "Y/W"                    ""            
## 
## $time
##  [1] 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958
## [15] 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972
## [29] 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986
## [43] 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000
## [57] 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010
## 
## $matrix
##     Yd T C H H_h Y N
## Yd   0 1 0 0   0 0 1
## T    0 0 0 0   0 0 1
## C    1 0 0 0   0 0 0
## H    0 1 0 0   0 0 0
## H_h  1 0 1 0   0 0 0
## Y    0 0 1 0   0 0 0
## N    0 0 0 0   0 1 0
## 
## $blocks
## $blocks[[1]]
## [1] 2 3 4 6 7 1 5
## 
## 
## attr(,"class")
## [1] "sfc"

Output data structure

The Gauss Seidel Algorithm

  1. Select initial values \(x^0\)
  2. While \(k<maxIter\) & \(\delta < tolValue\)
    1. For each \(i=1,...,n\): \[x_i^{k+1}=\frac{1}{a_{ii}}\left( b_i-\sum^{i-1}_{j=1}a_{ij} x_j^{k+1}-\sum_{j=i+1}^n a_{ij}x_j^k\right)\]
    2. Compute \(\delta\): \[\delta = \frac{x^{k+1}-x^k}{x^k}\]

System of (in)dependent equations

See Fenell et. al (2016)

Direct Acyclic graphs

simex<-sfc.model("Models/SIMEX.txt",modelName="SIMplest model with expectation")
layout(matrix(c(1,2),1,2))
plot.dag(sim,main="SIM" )
plot.dag(simex,main="SIMEX" )

In the case of a more complex model - Chapter 6

ch6 <- sfc.model("Models/ch6.txt",modelName="Chapter6_openmodel")
graphs = generate.DAG.collapse( adjacency = ch6$matrix )
par(mfrow = c(1,3))
# first plot the orgianl grpah
plot_graph_hierarchy( graphs$orginal_graph, main = "orginal graph" )
# plot hte nodes that form the strongly connected compoent
plot_graph_hierarchy( graphs$SCC_graph, main = "SCC nodes" )
# plot the result ing DAG when we take the condensation of the graph
plot_graph_hierarchy( graphs$DAG, main = "resulting DAG" )

Systems of dependent vs independent equations

#Doing the simulations
datasimex<-simulate(simex)
init = datasimex$baseline[66,]
simex<-sfc.addScenario(simex,"G",25,1960,2010,init)
datasimex<-simulate(simex)
datasim<-simulate(sim)
init = datasim$baseline[66,]
sim<-sfc.addScenario(sim,"G",25,1960,2010,init)
datasim<-simulate(sim)

#Plotting it all
layout(matrix(c(1,2),1,2))
matplot(sim$time,datasim$scenario_1[,c("H","C","Yd")],
    lty=1:3,type="l",xlab="",ylab="",main="SIM")
legend(x=1944,y=130,legend=c("Wealth","Consumption",
    "Disposable Income"), lty=1:3,bty="n")
matplot(simex$time,datasimex$scenario_1[,c("H","C","Yd",
    "Yd_e")],lty=1:4,type="l",xlab="",ylab="",main="SIMEX")
legend(x=1944,y=130,legend=c("Wealth","Consumption","Disposable Income",
   "Expecetd Disposable Income"),lty=1:4,bty="n")

Computational implications

Let’s see how much time it takes to run sim:

ptm <- proc.time()
data1<-simulate(sim)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],
                        "seconds"))
## [1] "Elapsed time is  3.67 seconds"

Now lets play with some of the parameters of the simulate function:

  1. tolValue
ptm <- proc.time()
data2<-simulate(sim,tolValue = 1e-3)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],"seconds"))
## [1] "Elapsed time is  0.176000000000002 seconds"
  1. maxIter
ptm <- proc.time()
data3<-simulate(sim, maxIter=10)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3],"seconds"))
## [1] "Elapsed time is  0.173999999999999 seconds"

Observing the results of the three simulations

kable(round(t(cbind(data1$baseline[c(1,2,20,40,66),c("Y")],
        data2$baseline[c(1,2,20,40,66),c("Y")],
        data3$baseline[c(1,2,20,40,66),c("Y")])),digits=3))
1945 1946 1964 1984 2010
NA 38.462 96.957 99.892 99.999
NA 38.464 96.874 99.577 99.911
NA 34.006 94.965 99.672 99.991

Block Gauss-Seidel

The order of equations matters, if first compute variables that do not depend on current period, this speeds the process. Define blocks of equation independent from the others.

print(simex$blocks)
## [[1]]
## [1] 7
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 6
## 
## [[4]]
## [1] 8
## 
## [[5]]
## [1] 2
## 
## [[6]]
## [1] 1 4
## 
## [[7]]
## [1] 5

Simulation of SIMEX

ptm <- proc.time()
dataex<-simulate(simex,tolValue = 1e-10)
print(paste("Elapsed time is ",proc.time()[3]-ptm[3]
                        ,"seconds"))
## [1] "Elapsed time is  0.117999999999999 seconds"

Results for SIMEX

kable(round(t(dataex$baseline[c(1,2,20,40,66),c("G","Y",
    "T","Yd","Yd_e","C","H","H_h")]),digits=3))
1945 1946 1964 1984 2010
G 20 20 20.000 20.000 20
Y NA 20 98.559 99.983 100
T NA 4 19.712 19.997 20
Yd 0 16 78.847 79.987 80
Yd_e NA 0 78.559 79.983 80
C NA 0 78.559 79.983 80
H 0 16 78.847 79.987 80
H_h 0 16 78.847 79.987 80

Output data structure

kable(dataex$baseline)
Yd T C H H_h Y Yd_e N alpha1 alpha2 theta G W iter block 1 iter block 2 iter block 3 iter block 4 iter block 5 iter block 6 iter block 7
1945 0.00000 NA NA 0.00000 0.00000 NA NA NA 0.6 0.4 0.2 20 1 0 0 0 0 0 0 0
1946 16.00000 4.00000 0.00000 16.00000 16.00000 20.00000 0.00000 20.00000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1947 28.80000 7.20000 16.00000 28.80000 28.80000 36.00000 16.00000 36.00000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1948 39.04000 9.76000 28.80000 39.04000 39.04000 48.80000 28.80000 48.80000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1949 47.23200 11.80800 39.04000 47.23200 47.23200 59.04000 39.04000 59.04000 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1950 53.78560 13.44640 47.23200 53.78560 53.78560 67.23200 47.23200 67.23200 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1951 59.02848 14.75712 53.78560 59.02848 59.02848 73.78560 53.78560 73.78560 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1952 63.22278 15.80570 59.02848 63.22278 63.22278 79.02848 59.02848 79.02848 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1953 66.57823 16.64456 63.22278 66.57823 66.57823 83.22278 63.22278 83.22278 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1954 69.26258 17.31565 66.57823 69.26258 69.26258 86.57823 66.57823 86.57823 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1955 71.41007 17.85252 69.26258 71.41007 71.41007 89.26258 69.26258 89.26258 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1956 73.12805 18.28201 71.41007 73.12805 73.12805 91.41007 71.41007 91.41007 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1957 74.50244 18.62561 73.12805 74.50244 74.50244 93.12805 73.12805 93.12805 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1958 75.60195 18.90049 74.50244 75.60195 75.60195 94.50244 74.50244 94.50244 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1959 76.48156 19.12039 75.60195 76.48156 76.48156 95.60195 75.60195 95.60195 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1960 77.18525 19.29631 76.48156 77.18525 77.18525 96.48156 76.48156 96.48156 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1961 77.74820 19.43705 77.18525 77.74820 77.74820 97.18525 77.18525 97.18525 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1962 78.19856 19.54964 77.74820 78.19856 78.19856 97.74820 77.74820 97.74820 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1963 78.55885 19.63971 78.19856 78.55885 78.55885 98.19856 78.19856 98.19856 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1964 78.84708 19.71177 78.55885 78.84708 78.84708 98.55885 78.55885 98.55885 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1965 79.07766 19.76942 78.84708 79.07766 79.07766 98.84708 78.84708 98.84708 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1966 79.26213 19.81553 79.07766 79.26213 79.26213 99.07766 79.07766 99.07766 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1967 79.40970 19.85243 79.26213 79.40970 79.40970 99.26213 79.26213 99.26213 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1968 79.52776 19.88194 79.40970 79.52776 79.52776 99.40970 79.40970 99.40970 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1969 79.62221 19.90555 79.52776 79.62221 79.62221 99.52776 79.52776 99.52776 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1970 79.69777 19.92444 79.62221 79.69777 79.69777 99.62221 79.62221 99.62221 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1971 79.75821 19.93955 79.69777 79.75821 79.75821 99.69777 79.69777 99.69777 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1972 79.80657 19.95164 79.75821 79.80657 79.80657 99.75821 79.75821 99.75821 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1973 79.84526 19.96131 79.80657 79.84526 79.84526 99.80657 79.80657 99.80657 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1974 79.87621 19.96905 79.84526 79.87621 79.87621 99.84526 79.84526 99.84526 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1975 79.90096 19.97524 79.87621 79.90096 79.90096 99.87621 79.87621 99.87621 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1976 79.92077 19.98019 79.90096 79.92077 79.92077 99.90096 79.90096 99.90096 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1977 79.93662 19.98415 79.92077 79.93662 79.93662 99.92077 79.92077 99.92077 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1978 79.94929 19.98732 79.93662 79.94929 79.94929 99.93662 79.93662 99.93662 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1979 79.95944 19.98986 79.94929 79.95944 79.95944 99.94929 79.94929 99.94929 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1980 79.96755 19.99189 79.95944 79.96755 79.96755 99.95944 79.95944 99.95944 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1981 79.97404 19.99351 79.96755 79.97404 79.97404 99.96755 79.96755 99.96755 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1982 79.97923 19.99481 79.97404 79.97923 79.97923 99.97404 79.97404 99.97404 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1983 79.98338 19.99585 79.97923 79.98338 79.98338 99.97923 79.97923 99.97923 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1984 79.98671 19.99668 79.98338 79.98671 79.98671 99.98338 79.98338 99.98338 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1985 79.98937 19.99734 79.98671 79.98937 79.98937 99.98671 79.98671 99.98671 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1986 79.99149 19.99787 79.98937 79.99149 79.99149 99.98937 79.98937 99.98937 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1987 79.99319 19.99830 79.99149 79.99319 79.99319 99.99149 79.99149 99.99149 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1988 79.99456 19.99864 79.99319 79.99456 79.99456 99.99319 79.99319 99.99319 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1989 79.99564 19.99891 79.99456 79.99564 79.99564 99.99456 79.99456 99.99456 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1990 79.99652 19.99913 79.99564 79.99652 79.99652 99.99564 79.99564 99.99564 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1991 79.99721 19.99930 79.99652 79.99721 79.99721 99.99652 79.99652 99.99652 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1992 79.99777 19.99944 79.99721 79.99777 79.99777 99.99721 79.99721 99.99721 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1993 79.99822 19.99955 79.99777 79.99822 79.99822 99.99777 79.99777 99.99777 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1994 79.99857 19.99964 79.99822 79.99857 79.99857 99.99822 79.99822 99.99822 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1995 79.99886 19.99971 79.99857 79.99886 79.99886 99.99857 79.99857 99.99857 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1996 79.99909 19.99977 79.99886 79.99909 79.99909 99.99886 79.99886 99.99886 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1997 79.99927 19.99982 79.99909 79.99927 79.99927 99.99909 79.99909 99.99909 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1998 79.99942 19.99985 79.99927 79.99942 79.99942 99.99927 79.99927 99.99927 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
1999 79.99953 19.99988 79.99942 79.99953 79.99953 99.99942 79.99942 99.99942 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2000 79.99963 19.99991 79.99953 79.99963 79.99963 99.99953 79.99953 99.99953 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2001 79.99970 19.99993 79.99963 79.99970 79.99970 99.99963 79.99963 99.99963 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2002 79.99976 19.99994 79.99970 79.99976 79.99976 99.99970 79.99970 99.99970 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2003 79.99981 19.99995 79.99976 79.99981 79.99981 99.99976 79.99976 99.99976 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2004 79.99985 19.99996 79.99981 79.99985 79.99985 99.99981 79.99981 99.99981 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2005 79.99988 19.99997 79.99985 79.99988 79.99988 99.99985 79.99985 99.99985 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2006 79.99990 19.99998 79.99988 79.99990 79.99990 99.99988 79.99988 99.99988 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2007 79.99992 19.99998 79.99990 79.99992 79.99992 99.99990 79.99990 99.99990 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2008 79.99994 19.99998 79.99992 79.99994 79.99994 99.99992 79.99992 99.99992 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2009 79.99995 19.99999 79.99994 79.99995 79.99995 99.99994 79.99994 99.99994 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2
2010 79.99996 19.99999 79.99995 79.99996 79.99996 99.99995 79.99995 99.99995 0.6 0.4 0.2 20 1 2 2 2 2 2 2 2

Checking the number of iteractions for SIM

kable(round(t(cbind(
    data1$baseline[c(1,2,20,40,66),c("iter block 1")],
    data2$baseline[c(1,2,20,40,66),c("iter block 1")],
    data3$baseline[c(1,2,20,40,66),c("iter block 1")]
    )),digits=3))
1945 1946 1964 1984 2010
0 293 218 175 119
0 89 15 2 1
0 10 10 10 10

Checking the number of iteractions for simex, no simulate parameters

kable(round(t(dataex$baseline[c(1,2,20,40,66),
    c("iter block 1","iter block 2","iter block 3",
    "iter block 4","iter block 5","iter block 6",
    "iter block 7")]),digits=3))
1945 1946 1964 1984 2010
iter block 1 0 2 2 2 2
iter block 2 0 2 2 2 2
iter block 3 0 2 2 2 2
iter block 4 0 2 2 2 2
iter block 5 0 2 2 2 2
iter block 6 0 2 2 2 2
iter block 7 0 2 2 2 2

Visualisation

Visualisation - Sankey Diagram

Flows in the UK Economy 2014, source: Burgess et al. (2016)

Flows in the UK Economy 2014, source: Burgess et al. (2016)

Visualisation - Graph representation

Example with models SIM-SIMEX

Graph representation of SIM and SIMEX, source: author’s computation

Graph representation of SIM and SIMEX, source: author’s computation

Shock propagation in more complex models

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Shock to CAR - Direct

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Shock to CAR - Lag 1: Interest rates

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Shock to CAR - Lag 2: Income and profits

Graph representation of Burgess et al. (2016)

Graph representation of Burgess et al. (2016)

Dynamical interaction with models

antoinegodin.shinyapps.io/SIMEX

Shiny application

Shiny application

Model with Portfolio choice

Balance Sheet

Households Production Government Central Bank Sum
Money +H -H 0
Bills +Bh -B +Bcb 0
Net worth -V +V 0
Sum 0 0 0 0 0

Transaction Flow Matrix

Households Production Government Central Bank Sum
Consumption -C +C 0
Gov. Exp. +G -G 0
Income = GDP +Y -Y 0
Interests +r(-1)*Bh(-1) -r(1)*B(-1) +r(-1)*Bcb(-1) 0
CB profits +r(-1)*Bcb(-1) -r(1)*Bcb(-1) 0
Taxes -T +T 0
Change in Money -\(\Delta\) H +\(\Delta\) H 0
Change in Bills -\(\Delta\) Bh +\(\Delta\) B -\(\Delta\) Bcb 0
Sum 0 0 0 0 0

Behavioural Equations (from PC.txt file)

#######Determination of output - eq. 4.1
y = cons + g
#######Disposable income - eq. 4.2
yd = y - t + r(-1)*b_h(-1)
#######Tax payments - eq. 4.3
t = theta*(y + r(-1)*b_h(-1))
#######Wealth accumulation - eq. 4.4
v = v(-1) + (yd - cons)
#######Consumption function - eq. 4.5
cons = alpha1*yd + alpha2*v(-1)
#######Cash money - eq. 4.6
h_h = v - b_h  
#######Demand for government bills - eq. 4.7
b_h = v*(lambda0 + lambda1*r - lambda2*(yd/v))
#######Supply of government bills - eq. 4.8
b_s = b_s(-1) + (g + r(-1)*b_s(-1)) - (t + r(-1)*b_cb(-1))
#######Supply of cash - eq. 4.9
h_s = h_s(-1) + b_cb - b_cb(-1)
#######Government bills held by the central bank - eq. 4.10
b_cb = b_s - b_h
#######Interest rate as policy instrument - eq. 4.11
r = r_bar

Portfolio Equations

How to allocate wealth between bills and money?

  1. Quantity theory of money: links money balances to the flow of income
  2. Liquidity preferences: money balances some proportion of total wealth

Merging the two into Tobin’s portfolio equations:

\[\begin{align} \frac{H_h}{V}&=(1-\lambda_0)-\lambda_1\cdot r+\lambda_2\cdot\left(\frac{YD}{V}\right)\tag{4.6A}\\ \frac{B_h}{V}&=\lambda_0+\lambda_1\cdot r-\lambda_2\cdot\left(\frac{YD}{V}\right)\tag{4.7}\\ H_h&=V-B_h\tag{4.6} \end{align}\]

Setting up the environment

Before doing any modelling, we need to load the package in the R environment.

library(PKSFC)

Then, you need to download the attached ‘PC.txt’ file and save it in the folder of your choice. Make sure to set the working directory where you saved the downloaded file. In command line this looks like this but if you use Rstudio, you can use the graphical interface as well (Session>Set Working Directory>Choose Directory)

setwd("pathToYourDirectory")

Loading the model

The first thing to do is to load the model and check for completeness.

pc<-sfc.model("Models/PC.txt",modelName="Portfolio Choice Model")
pc<-sfc.check(pc,fill=FALSE)

Let’s have a look at the graphical representation of model PC. In order to do so, we need to load the Rgraphviz library:

library("Rgraphviz")

We can now look at the graph of model PC:

plot.dag(pc,main="PC" )

You can see that there is a cycle in the graph, implying that GDP, taxes, disposable income and consumption are determined all together (and that they fully adapt to any shock applied to the economy). While this is a mathematical property of the system of equations representing the economy we wish to model, it has an economic meaning and you want to be sure that this is what you believe to be the best representation of what you have in mind.

We are now ready to simulate the model

datapc<-simulate(pc)

Expectations and random shocks

The first of experiment is meant to observe the buffer role of money, in case of random shocks applied to expected disposable income. To do these, we need to modify slightly model pc and change the equations determining consumption, demand for bonds and money, and the expectations on income and wealth. We will call the new model pcRand.

pcRand<-sfc.editEqus(pc,list(list(var="cons",eq="alpha1*yde+alpha2*v(-1)"),
      list(var="b_h",eq="ve*(lambda0 + lambda1*r - lambda2*(yde/ve))")))
pcRand<-sfc.addEqus(pcRand,list(
  list(var="yde",eq="yd(-1)*(1+rnorm(1,sd=0.1))",
         desc="Expected disposable income depending on random shocks"),
  list(var="h_d",eq="ve-b_h",desc="Money demand"),
  list(var="ve",eq="v(-1)+yde-cons",desc="Expected disposable income")))
pcRand<-sfc.editVar(pcRand,var="yd",init=86.48648)
pcRand<-sfc.check(pcRand)

You probably have noticed that the yde equation contains the R function rnorm which allows you to extract a random number from a normal distribution centred around 0, with standard variation of 0.1. The package indeed allows you to use any R function such as min, max or rnorm. However, because of the way the Gauss-Seidel algorithm works, we are facing a problem. Indeed, as the rnorm function extract a number in each iteration, this will mean that the algorithm cannot converge since the difference between each iteration of the algorithm depends of the difference between each extraction. This is why we need to restrict the number of iteration of the Gauss-Seidel.

Before simulating the model, let’s have a look at how the graph of the model has changed:

plot.dag(pcRand,main="PC Random" )

You can now see that the cycle observed in the original model has disappeared. Indeed, consumption doesn’t depend on disposable income any more but on expected disposable income. Let’s now simulate the model.

datapcRand<-simulate(pcRand,maxIter=2)

This replicates figure 4.1, page 110

plot(pcRand$time,datapcRand$baseline[,"h_h"],type="l",xlab="",ylab="",lty=1,lwd=2,
     ylim=c(1*min(datapcRand$baseline[,c("h_h","h_d")],na.rm=T),
            1.2*max(datapcRand$baseline[,c("h_h","h_d")],na.rm=T)))
lines(pcRand$time,datapcRand$baseline[,"h_d"],lty=2,lwd=2)
legend(x=1950,y=1.2*max(datapcRand$baseline[,c("h_h","h_d")],na.rm=T),
       legend=c("Money held","Money demand"),lty=c(2,1),lwd=2,bty="n")

This graph highlights the buffer role of certain stocks in PK-SFC models. Indeed, because expectations are incorrect or because demand might not be equal to supply in any market, at least one stock will not be equal to the targeted level. As highlighted by Foley (1975), in a model without perfect foresight you need a buffer stock in order to obtain equilibrium between demand and supply. The role of buffer stocks in PK-SFC model is thus fundamental and is at the hart of the approach used by Godley (1999) in his seven unsustainable processes. It is by observing the dynamics of certain stock-flow norms that you are able to observe the unsustainable processes evolving in an economy, because stocks precisely absorb disequilibrium.

Interest rates impacts

The graphs presented on the paper are based on the model PCEX, which includes expectation on disposable income and wealth. We will first create the model before simulating it.The shock represents an 100 basis point increase in interest rates.

pcex<-sfc.editEqus(pc,list(
  list(var="cons",eq="alpha1*yde+alpha2*v(-1)"),
  list(var="b_h",eq="ve*(lambda0 + lambda1*r - lambda2*(yde/ve))")))
pcex<-sfc.addEqus(pcex,list(
  list(var="yde",eq="yd(-1)",desc="Expected disposable income"),
  list(var="h_d",eq="ve-b_h",desc="Money demand"),
  list(var="ve",eq="v(-1)+yde-cons",desc="Expected disposable income")))
pcex<-sfc.editVar(pcex,var="yd",init=86.48648)
pcex<-sfc.check(pcex)
init = datapc$baseline[56,]
pcex<-sfc.addScenario(model=pcex,vars="r_bar",values=0.035,inits=1960,
                                            ends=2010,starts=init)
datapcex<-simulate(pcex)

This replicates plot figure 4.3. p. 112

time2=c(1950:2000)
plot(time2, datapcex$scenario[as.character(time2),"h_h"]/
            datapcex$scenario[as.character(time2),"v"],
         type="l",xlab="",ylab="",lty=1,lwd=2)
par(new=T)
plot(time2,datapcex$scenario[as.character(time2),"b_h"]/
            datapcex$scenario[as.character(time2),"v"],
     lty=2, lwd=2,type="l",axes=F,ylab="",xlab="")
axis(4,pretty(c(0.750, 1.1*max(datapcex$scenario[,"b_h"]/
            datapcex$scenario[,"v"],na.rm=T))))
legend(x=1970,y=0.78,legend=c("Share of Bills",
        "Share of Money"),lty=c(2,1),lwd=2,bty="n")
grid()

No surprise here, the change in return rates coming from the shock led to a re-allocation between Bills and Money, in favour of the former.

This replicates fig 4.4 . p 113

time2=c(1950:2000)
matplot(time2,datapcex$scenario[as.character(time2),c("yd","cons")]
                ,type="l",xlab="",ylab="",lty=1:2,col=1,lwd=2)
legend(x=1970,y=88,legend=c("Disposable Income","Consumption"),
             lwd=c(2,2),lty=c(1,2),bty="n")
grid()

These results are more surprising as increased interest rates lead to an increase both in disposable income and consumption. To understand this, one has to bear in mind that increase in interest rates leads to an increase in government spending and thus to an increase in the fiscal stance which determines the steady state:

\(Y^\star=\frac{G+r \cdot B^{\star}_h \cdot (1-\theta)}{\theta}\)

Endogenous propensities to consume

We change the propensities to consume to incorporate the fact that they might be impacted by interest rates and thus have

\[\begin{equation} \alpha_1=\alpha_10-\iota\cdot r_{-1} \tag{4.32} \end{equation}\]

This piece of code shows how to implement this

pc_end<-sfc.addEqus(pcex,list(
  list(var="alpha1",eq="alpha10-iota*r(-1)")))
pc_end<-sfc.addVars(pc_end,list(
  list(var="alpha10",init="0.7",
         desc="Endoegnous propensity to consume - autonomous term"),
  list(var="iota",init="10",
         desc="Endoegnous propensity to consume - interest rate impact")))
pc_end$scenarios<-NULL
pc_end<-sfc.check(pc_end)
datapc_end<-simulate(pc_end)
init = datapc_end$baseline[56,]
pc_end<-sfc.addScenario(model=pc_end,vars="r_bar",values=0.035,
                                                inits=1960,ends=2010,starts=init)
datapc_end<-simulate(pc_end)

This replicates plot figure 4.9. p. 123

time2=c(1950:2000)
matplot(time2,
     datapc_end$scenario_1[as.character(time2),c("v","yd","cons","y")],
     type="l",xlab="",ylab="",lty=1:4,lwd=2,col=1)
grid()
legend(x=1950,y=100,legend=c("Households wealth","Disposable income",
            "Consumption","GDP"),lty=1:4,lwd=2,bty="n")

This shows that now the short-term dynamics display the Keynesian paradox of thrift but the long-run steady state still shows a positive impact of interests rates on GDP. The short-run recession is governed by the decrease in the propensity to consumed and is slowly compensated by the increase in wealth (and the increased public debt as a counterpart) leading to the new steady state. The steady state is still driven by the fiscal stance.

This replicates plot figure 4.10. p. 123

timelag=c(1950:2000)
time2=c(1951:2001)
plot(time2,
     datapc_end$scenario_1[as.character(time2),c("g")]+
            datapc_end$scenario_1[as.character(timelag),c("r")]*
            datapc_end$scenario_1[as.character(timelag),c("b_s")]-
            datapc_end$scenario_1[as.character(timelag),c("r")]*
            datapc_end$scenario_1[as.character(timelag),c("b_cb")],
     type="l",xlab="",ylab="",lty=1:4,lwd=2,col=1,ylim=c(20.8,23))
lines(time2,datapc_end$scenario_1[as.character(time2),"t"],type="l",xlab="",
            ylab="",lty=2,lwd=2)
grid()
legend("topleft",legend=c("Government expenditures plus net debt services",
            "Tax revenues"),lty=1:4,lwd=2,bty="n")

Debt to GDP at the steady state

\[\begin{equation} \frac{V^\star}{Y^\star}=\frac{\frac{1-\alpha_1}{\alpha_2}}{1+\left[\frac{\theta}{1-\theta}\right]-r\cdot \left[(\lambda_0+\lambda_1\cdot r)\cdot \frac{1-\alpha_1}{\alpha_2}-\lambda_2\right]}\tag{4.33} \end{equation}\]

The Maastricht treaty: The reference values referred to […] are:

Liquidity Preferences and model LP

Value of perpetuity, interest rates, exepcted returns and capital gains

The value of a financial asset is supposed to reflect the net present value of future cash flows. Thus for a long-term bond who is never redeemed, the price \(p_{bL}\) of the bond is given by

\(p_{bL} = \sum \frac{1}{(1+r_{bL})^t} = \frac{1}{r_{bL}}\)

Return rate of an asset is equal to both the yield (interests, dividends) and the capital gains normalized to the nominal value of the asset. In the case of a long-term bond:

\(Rr_{bL} = r_{bL}(-1)+\frac{\Delta p_{bL}}{p_{bL}(-1)} = \frac{1+p_{bL}-p_{bL}(-1)}{p_{bL}(-1)}\)

However, what matter for the households when making their decision is the expected price of bonds \(p_{bL}^e\), given the current price of bonds. This leads to the pure expected rate of return.

\(PERr_{bL} = r_{bL} + \frac{p^e_{bL}-p_{bL}}{p_{bL}}\)

Because expectation are not followed by everyone, households might incorporate a weight \(\xi\) into their expectation formation. This leads to the expected return rate

\(ERr_{bL} = r_{bL} + \xi \frac{p^e_{bL}-p_{bL}}{p_{bL}}\)

The change in prices of bonds leads to capital gains which is one of two parts of the change in nominal value of financial assets. More precisely:

\[\begin{align*} \Delta (p_{BL}\cdot BL) &=\left(p_{bL} \cdot BL\right) - \left(p_{bL}(-1) \cdot BL(-1)\right)\\ &= p_{BL} \cdot BL-p_{BL}\cdot BL(-1) + p_{BL}\cdot BL(-1)-p_{bL}(-1) \cdot BL(-1)\\ &= p_{BL} \left( BL-BL(-1)\right) + BL(-1) \left(p_{BL}-p_{BL}(-1)\right)\\ &= p_{BL} \cdot \Delta BL+ BL(-1)\cdot \Delta p_{BL} \end{align*}\]

Thus total changes in nominal values of an asset is equal to the quantity of assets bought (or sold) at the current price (\(p_{BL} \cdot \Delta BL\)) and to the capital gains (or losses) made on the stock of assets held at the beginning of the period (\(BL(-1)\cdot \Delta p_{BL}\))

Ostergaard diagram

Ostergaard diagram

Balance Sheet

Households Production Government Central Bank Sum
Money +H -H 0
Bills +Bh -B +Bcb 0
Bonds +BL.pbL -BL.pbL 0
Net worth -V +V 0
Sum 0 0 0 0 0

Transaction Flow Matrix

Households Production Government Central Bank Sum
Consumption -C +C 0
Gov. Exp. +G -G 0
Income = GDP +Y -Y 0
Interests on bills +rb(-1)*Bh(-1) -rb(1)*B(-1) +rb(-1)*Bcb(-1) 0
Interests on bonds +BL(-1) -BL(-1) 0
CB profits +r(-1)*Bcb(-1) -r(1)*Bcb(-1) 0
Taxes -T +T 0
Change in Money -\(\Delta\) H +\(\Delta\) H 0
Change in Bills -\(\Delta\) Bh +\(\Delta\) B -\(\Delta\) Bcb 0
Change in Bonds -\(\Delta\) BL.pbL +\(\Delta\) BL.pbL 0
Sum 0 0 0 0 0
Memo: Capital gains -\(\Delta\) pbL.BL(-1) +\(\Delta\) pbL.BL(-1) 0

Equation list

# MODEL
# Determination of output - eq. 5.1
y = cons + g
# Regular disposable income - eq. 5.2
yd_r = y - t + r_b(-1)*b_h(-1) + bl_h(-1)
# Tax payments - eq. 5.3
t = theta*(y + r_b(-1)*b_h(-1) + bl_h(-1))
# Wealth accumulation - eq. 5.4
v = v(-1) + (yd_r - cons) + cg
# Capital gains on bonds - eq. 5.5
cg = (p_bl-p_bl(-1))*bl_h(-1)
# Consumption function - eq. 5.6
cons = alpha1*yd_r_e + alpha2*v(-1)
# Expected wealth - eq. 5.7
v_e = v(-1) + (yd_r_e - cons) + cg
# Cash money - eq. 5.8
h_h = v - b_h - p_bl*bl_h
# Demand for cash - eq. 5.9
h_d = v_e - b_d - p_bl*bl_d
# Demand for government bills - eq. 5.10
b_d = v_e*(lambda20 + lambda22*r_b - lambda23*er_rbl - lambda24*(yd_r_e/v_e))
# Demand for government bonds - eq. 5.11
bl_d = v_e*(lambda30 - lambda32*r_b + lambda33*er_rbl - lambda34*(yd_r_e/v_e))/p_bl
# Bills held by households - eq. 5.12
b_h = b_d
# Bonds held by households - eq. 5.13
bl_h = bl_d
# Supply of government bills - eq. 5.14
b_s = b_s(-1) + (g + r_b(-1)*b_s(-1) + bl_s(-1)) - (t + r_b(-1)*b_cb(-1)) -p_bl*(bl_s-bl_s(-1))
# Supply of cash - eq. 5.15
h_s = h_s(-1) + b_cb - b_cb(-1)
# Government bills held by the central bank - eq. 5.16
b_cb = b_s - b_h
# Supply of government bonds - eq. 5.17
bl_s = bl_h
# Expected rate of return on bonds - eq. 5.18
er_rbl = r_bl+chi*(p_bl_e - p_bl)/p_bl
# Interest rate on bonds - eq. 5.19
r_bl = 1/p_bl
# Expected price of bonds - eq. 5.20
p_bl_e = p_bl
# Expected capital gains - eq. 5.21
cg_e = chi*(p_bl_e - p_bl)*bl_h
# Expected regular disposable income - eq. 5.22
yd_r_e = yd_r(-1)
# Interest rate on bills - eq. 5.23
r_b = r_bar
# Price of bonds - eq. 5.24
p_bl = p_bl_bar

Let’s have a look at the graph of model LP

lp<-sfc.model(fileName="Models/LP.txt",modelName="Model Liquidity Preference")
plot.dag(lp)

The adding up of the portfolio behavior

The vertical conditions state that the sum of the exogenous components sums to 1, i.e. you cannot allocate more than the wealth and that the effect of the various return rates and disposable income are subtracting the holding of one asset to allow the increase of the holding of another.

\[\begin{align} \lambda_{10} + \lambda_{20} + \lambda_{30} &= 1 \tag{ADUP.1}\\ \lambda_{11} + \lambda_{21} + \lambda_{31} &= 0 \tag{ADUP.2}\\ \lambda_{12} + \lambda_{22} + \lambda_{32} &= 0 \tag{ADUP.3}\\ \lambda_{13} + \lambda_{23} + \lambda_{33} &= 0 \tag{ADUP.4}\\ \lambda_{14} + \lambda_{24} + \lambda_{34} &= 0 \tag{ADUP.5} \end{align}\]

The horizontal conditions, added by Godley (1996) state that the impact of an increase in the ‘own rate’ of an asset should be equal to the impact of a fall in all the other rates.

\[\begin{align} \lambda_{11} &= - \left( \lambda_{12} + \lambda_{13}\right) \tag{ADUP.6}\\ \lambda_{22} &= - \left( \lambda_{21} + \lambda_{23}\right) \tag{ADUP.7}\\ \lambda_{33} &= - \left( \lambda_{31} + \lambda_{32}\right) \tag{ADUP.8} \end{align}\]

Other authors propose the symmetry conditions that state that an increase in the return rate of one asset A will have the same impact on the holding of asset B as the increase in the return rate of asset B will have on the holding of asset A

\[\begin{align} \lambda_{12} &= \lambda_{21} \tag{ADUP.9}\\ \lambda_{13} &= \lambda_{31} \tag{ADUP.10}\\ \lambda_{23} &= \lambda_{32} \tag{ADUP.11} \end{align}\]

Note that vertical + symmetry imply automatically horizontal but that vertical + horizontal does not necessarily imply symmetry.

Higher interest rates

The first scenario run looks at the impact of interest rates (short- and long-term) on real demand. We assume that the government increases short-term rates from 3% to 4% and the long-term rate from 5% to 6.66%. We assume that this increase is completely unexpected and that the Treasury is believed in that no other change will occur.

lp<-sfc.addScenario(model=lp,vars=list(c("r_bar","p_bl_bar")),values=list(c(0.04,15)),
                            inits=1960,ends=2000)
datalp<-simulate(lp)

This replicates plot figure 5.2. p. 152

time2=c(1950:2000)
plot(time2, datalp$scenario_1[as.character(time2),"v"]/
            datalp$scenario[as.character(time2),"yd_r"],
         type="l",xlab="",ylab="",lty=1,lwd=2)
legend("bottomright",legend=c("Wealth to income ratio"),lty=c(1),lwd=2,bty="n")
grid()

We see that the increase in long-term interest rate leads to capital losses on bonds, leading to a decrease in wealth to income ratio. However because interest rates have increase, households have now a higher disposable income leading to a replenishing of their wealth. Note that as in model PC, the steady state GDP is given by

\(Y^\star=\frac{G+\left(r_B\cdot B^{\star}_h+BL_S^\star\right) \cdot (1-\theta)}{\theta} = \frac{G_{NT}}{\theta}\)

This replicates plot figure 5.3. p. 152

time2=c(1950:2000)
matplot(time2, datalp$scenario_1[as.character(time2),c("yd_r","cons")],
         type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("topleft",legend=c("Disposable income","Consumption"),lty=c(1,2),lwd=2,bty="n")
grid()

This replicates plot figure 5.4. p. 153

time2=c(1950:2000)
matplot(time2, cbind(datalp$scenario_1[as.character(time2),c("b_h")],
            datalp$scenario_1[as.character(time2),c("p_bl")]*
            datalp$scenario_1[as.character(time2),c("bl_h")])/
                datalp$scenario_1[as.character(time2),c("v")],
         type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("center",legend=c("Bonds to wealth ratio","Bills to wealth ratio"),
        lty=c(2,1),lwd=2,bty="n")
grid()

Introducing household liquidity

In order to introduce liquidity preferences, we need to modify the model. The first modification is to change the expectation on bonds price to include (i) a learning element and (ii) a stochastic term. Because expectation are going to play a role, we also need to modify the price equation for bonds such that it now reflects a corridor approach to bond pricing. The assumption is thus that the Treasury will try to pin down the price of bonds but only to a certain extent. If the demand is too large (or too low) the treasury will let prices (and interest) float. We assume that the government aims at a certain debt structure in terms of maturity.

#Adding the equations
lp2<-sfc.addEqus(lp,list(
    list(var="z1",equ = "TP>top"),
    list(var="z2",equ = "TP<bot"),
    list(var="TP",equ = "(bl_h(-1)*p_bl(-1))/(bl_h(-1)*p_bl(-1)+b_h(-1))")))
#Adding the parameters
lp2<-sfc.addVars(lp2,list(
    list(var="betae",init=0.5),
    list(var="beta",init=0.01),
    list(var="add",init=0),
    list(var="top",init=0.495),
    list(var="bot",init=0.505)
))
#Adding initial values to the endogenous variables that have lags
lp2<-sfc.editEnd(lp2,var="p_bl_e", init=20)
lp2<-sfc.editEnd(lp2,var="p_bl", init=20)
#Modifying the existing equations
lp2<-sfc.editEqus(lp2,list(
    list(var="p_bl", eq="(1+z1*beta-z2*beta)*p_bl(-1)"),
    list(var="p_bl_e",equ = "p_bl_e(-1)-betae*(p_bl_e(-1) - p_bl) + add")))
#Removing all existing scenarios (coming from LP)
lp2$scenarios<-NULL
#Checking that the model is complete
lp2<-sfc.check(lp2,fill=F)
#Adding the scenario
lp2<-sfc.addScenario(model=lp2,vars=list(c("r_bar")),values=list(c(0.035)),inits=1960,
                                            ends=2000)
lp2<-sfc.addScenario(model=lp2,vars=list(c("add")),values=list(c(-3)),inits=c(1955),
                                            ends=1955)

Let’s have a look at the graph of model LP2

plot.dag(lp2)

We are now ready to simulate the model

datalp2<-simulate(lp2)

This replicates plot figure 5.5. p. 156

time2=c(1950:2000)
plot(time2,datalp2$scenario_1[as.character(time2),"r_b"],
     lty=2, lwd=2,type="l",ylab="",xlab="",ylim=c(0.029,0.036))
par(new=T)
plot(time2, datalp2$scenario_1[as.character(time2),"r_bl"],
         type="l",xlab="",axes=F,ylab="",lty=1,lwd=2,ylim=c(0.0495,0.052))
axis(4,pretty(c(0.0495,0.052)))
legend("bottomright",legend=c("Long term interest rate",
        "Short term interest rate"),lty=c(1,2),lwd=2,bty="n")
grid()

This replicates plot figure 5.6. p. 156

time2=c(1950:2000)
matplot(time2, datalp2$scenario_1[as.character(time2),c("TP","bot","top")],
         type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend(x=1950,y=0.502,legend=c("Share of bonds in government debt held by households"),
             lty=c(2),lwd=2,bty="n")
grid()

Let’s add a second scenario by which there is an expected fall in the price of long-term bonds.

lp2<-sfc.addScenario(model=lp2,vars=list(c("add")),values=list(c(-3)),inits=c(1955),
                                            ends=1955)
datalp2<-simulate(lp2)

This replicates plot figure 5.7. p. 157

time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("r_bl")],
         type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend(x=1950,y=0.502,legend=c("Share of bonds in government debt held by households"),
             lty=c(2),lwd=2,bty="n")
grid()

This replicates plot figure 5.8. p. 158

time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("p_bl_e","p_bl")],
         type="l",xlab="",ylab="",lty=c(1,2), col=1,lwd=2)
legend("center",legend=c("Expected price of bonds","Actual price of bonds"),
             lty=c(1,2),lwd=2,bty="n")
grid()

This replicates plot figure 5.9. p. 158

time2=c(1950:2000)
matplot(time2, datalp2$scenario_2[as.character(time2),c("TP","bot","top")],
         type="l",xlab="",ylab="",lty=c(2,1,1), col=1,lwd=2)
legend("center",legend=c("Share of bonds in government debt held by households"),
             lty=c(2),lwd=2,bty="n")
grid()

Modelling the firm: price, investment and profits

The role of banks: credit creation and investment

Functional income distribution

Consumption function

Steady state of BMW

#Supply of consumption goods - eq. 7.1
 c_s = c_d
#Supply of investment goods - eq. 7.2
 i_s = i_d
#Supply of labour - eq. 7.3
 n_s = n_d
#Transactions of the firms
#GDP - eq. 7.5
 y = c_s + i_s
#Wage bill - eq. 7.6
 wb_d = y - r_l(-1)*l_d(-1) - af
#Depreciation allowances - eq. 7.7
 af = delta*k(-1)
#Demand for bank loans - eq. 7.8
 l_d = l_d(-1) + i_d - af
#Transactions of households
#Disposable income - eq. 7.9
 yd = wb_s + r_m(-1)*m_h(-1)
#Bank deposits held by households - eq. 7.10
 m_h = m_h(-1) + yd - c_d
#The wage bill
#"Supply" of wages - eq. 7.13
 wb_s = w*n_s
#Labour demand - eq. 7.14
 n_d = y/pr
#Wage rate - eq. 7.15
 w = wb_d/n_d
#Household behaviour
#Demand for consumption goods - eq. 7.16
 c_d = alpha0 + alpha1*yd + alpha2*m_h(-1)
#The investment behaviour
#Accumulation of capital - eq. 7.17
 k = k(-1) + i_d - da
#Depreciation allowances - eq. 7.18
 da = delta*k(-1)
#Capital stock target - eq. 7.19
 k_t = kappa*y(-1)
# below is an additional equation I've defined to get output to capital ratio as in figure 7.4, the model will run even without equation
OCR = y/k(-1)  
#Demand for investment goods - eq. 7.20
 i_d = gamma*(k_t - k(-1)) + da

Out of equilibrium values

BMW<-sfc.model("Models/BMW.txt")
## Warning in sfc.check(model, fill = fill): The following variables have lags
## but no initial values: - r_l - r_m
plot.dag(BMW)

Stability

More on profits and inventories

Components Firms Current account Firms Capital account
Sales \(+S\)
Change in the value of inventories +\(\Delta IN\) -\(\Delta IN\)
Wages \(-WB\)
Interest on loans \(-r_{l-1} L_{-1}\)
Entrepreneurial profits \(-F\)
Change in loans \(+\Delta L\)
Sum 0 0

References

Foley, D. 1975. “On Two Specifications of Asset Equilibrium in Macroeconomic Models.” Journal of Political Economy 83-2 (2): 303–24.

Godley, Wynne. 1999. “Seven Unsustainable Processes: Medium-Term Prospects and Policies for the United States and the World.” Strategic Analysis. The Levy Economic Institute of Bard College.